Note: install.packages(c(“tidycensus”, “tidyverse”, “sf”, “mapview”, “geepack”))
Note: You only need to do this once per device using install = TRUE.
##
## Call:
## glm(formula = N_SHOOTINGS_2025 ~ N_BROKEN_LIGHTS_2025, family = poisson,
## data = final_cross_section_2025)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.91705 0.19716 -9.723 < 2e-16 ***
## N_BROKEN_LIGHTS_2025 0.04535 0.01503 3.018 0.00254 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 252.75 on 309 degrees of freedom
## Residual deviance: 244.93 on 308 degrees of freedom
## AIC: 365.33
##
## Number of Fisher Scoring iterations: 6
## (Intercept) N_BROKEN_LIGHTS_2025
## 0.1470408 1.0463913
##
## Call:
## geepack::geeglm(formula = N_SHOOTINGS ~ N_BROKEN_LIGHTS + created_year,
## family = poisson, data = final_longitudinal_data, id = GEOID,
## corstr = "unstructured")
##
## Coefficients:
## Estimate Std.err Wald Pr(>|W|)
## (Intercept) 448.829995 56.799488 62.442 2.78e-15 ***
## N_BROKEN_LIGHTS 0.012353 0.004656 7.039 0.00798 **
## created_year -0.222397 0.028100 62.639 2.44e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation structure = unstructured
## Estimated Scale Parameters:
##
## Estimate Std.err
## (Intercept) 2.459 0.2596
## Link = identity
##
## Estimated Correlation Parameters:
## Estimate Std.err
## alpha.1:2 0.7595 0.09163
## alpha.1:3 0.4816 0.08234
## alpha.1:4 0.4816 0.07975
## alpha.1:5 0.5573 0.08891
## alpha.1:6 0.2690 0.06500
## alpha.2:3 0.6232 0.09262
## alpha.2:4 0.6866 0.08908
## alpha.2:5 0.6844 0.10104
## alpha.2:6 0.3365 0.08126
## alpha.3:4 0.4244 0.06458
## alpha.3:5 0.5023 0.09885
## alpha.3:6 0.2154 0.05049
## alpha.4:5 0.4427 0.05812
## alpha.4:6 0.2284 0.06247
## alpha.5:6 0.2165 0.05758
## Number of clusters: 310 Maximum cluster size: 6
## (Intercept) N_BROKEN_LIGHTS created_year
## 8.402e+194 1.012e+00 8.006e-01
final_cross_section_2025_DAY <- manhattan_pop |>
left_join(shootings_per_tract_25_DAY, by = "GEOID") |>
left_join(lights_per_tract_25, by = "GEOID") |> # lights_per_tract_25 calculated in Part 2
mutate(
N_BROKEN_LIGHTS_2025 = replace_na(N_BROKEN_LIGHTS_2025, 0),
N_SHOOTINGS_2025 = replace_na(N_SHOOTINGS_2025, 0)
) |>
st_drop_geometry()
poisson_model_2025_DAY <- glm(
N_SHOOTINGS_2025 ~ N_BROKEN_LIGHTS_2025,
family = poisson,
data = final_cross_section_2025_DAY
)
summary(poisson_model_2025_DAY)
print(exp(coef(poisson_model_2025_DAY)))
final_cross_section_2025_NIGHT <- manhattan_pop |>
left_join(shootings_per_tract_25_NIGHT, by = "GEOID") |>
left_join(lights_per_tract_25, by = "GEOID") |>
mutate(
N_BROKEN_LIGHTS_2025 = replace_na(N_BROKEN_LIGHTS_2025, 0),
N_SHOOTINGS_2025 = replace_na(N_SHOOTINGS_2025, 0)
) |>
st_drop_geometry()
poisson_model_2025_NIGHT <- glm(
N_SHOOTINGS_2025 ~ N_BROKEN_LIGHTS_2025,
family = poisson,
data = final_cross_section_2025_NIGHT
)
summary(poisson_model_2025_NIGHT)
print(exp(coef(poisson_model_2025_NIGHT)))
| Time Period | IRR (Incidence Rate Ratios) | P-value |
|---|---|---|
| Daytime (6 AM - 7:59 PM) | 1.057 | 0.0081 |
| Nighttime (8 PM - 5:59 AM) | 1.036 | 0.0973 |
poisson_gee_model_DAY <- geepack::geeglm(
N_SHOOTINGS ~ N_BROKEN_LIGHTS + created_year,
data = final_longitudinal_data_DAY,
id = GEOID,
family = poisson,
corstr = "unstructured"
)
summary(poisson_gee_model_DAY)
print(exp(coef(poisson_gee_model_DAY)))
poisson_gee_model_NIGHT <- geepack::geeglm(
N_SHOOTINGS ~ N_BROKEN_LIGHTS + created_year,
data = final_longitudinal_data_NIGHT,
id = GEOID,
family = poisson,
corstr = "unstructured"
)
summary(poisson_gee_model_NIGHT)
print(exp(coef(poisson_gee_model_NIGHT)))
| Time Period | IRR (Incidence Rate Ratios) | P-value (Wald Test) |
|---|---|---|
| Daytime (6 AM - 7:59 PM) | 1.01 | 0.058 |
| Nighttime (8 PM - 5:59 AM) | 1.02 | 0.000 |